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Abstract 

Optimal mass transport is described by an approximation of transport cost via semi-discrete 
costs. The notions of optimal partition and optimal strong partition are given as well. We also 
suggest an algorithm for computation of Optimal Transport for general cost functions induced 
by an action, an asymptotic error estimate and several numerical examples of optimal partitions. 


1 Introduction 


Optimal mass transport (OMT) goes back to the pioneering paper of Monge 15 at the 18th century. 
In 1942, L. Kantorovich [T^ observed that OMT can be relaxed into an infinite dimensional linear 


programming in measure spaces. As such , it has a dual formulation which is very powerful and was 
later (1987) used by Brenier to develop the theory of Polar factorization of positive measures. 
OMT has many connections with PDE, kinetic theory, fluid dynamics, geometric inequalities, 
probability and many other fields in mathematics as well as in computer science and economy. 


Even though finite dimensional (or discrete) OMT is well understood, its extension to infinite 
dimensional measure spaces poses a great challenge, e.g. uniqueness and regularity theory of fully 
non-linear PDE such as the Monge-Amper equation [^. 

We suggest to investigate a bridge between finite (’’discrete”) and infinite (’’continuum”) di¬ 
mensional OMT. This notion of semi-discrete OMT leads naturally to optimal partition of measure 
spaces. Our motivation in this paper is the development of numerical method for solving OMT. 
Efficient algorithms are of great interest to many fields in operational research and, recently, also 
for optical design and computer vision (” earth moving metric”) . 

When dealing with numerical approximations for OMT, the problem must be reduced to a 
discrete, finite OMT (with, perhaps, very large number of degrees of freedom). Discrete OMT is 
often called the assignment problem. This is, in fact, a general title for a variety of linear and 
quadratic programming. It seems that the first efficient algorithm was the so called ’’Hungarian 
Algorithm”, after two Hungarian mathematicians. See and the survey paper 

for many other relevant references. 


The deterministic, finite assignment problem is easy to formulate. We are given n men and n 
women. The cost of matching man z to a woman j is Cjj. The object is to find the assignment 
(matching) i ^ j, given in terms of a permutation j = T(i) which minimize the total cost of 
matching 

When replacing the deterministic assignment by a probabilistic one, we assign the probability 
pi > 0 for matching man i to woman j. The discrete assignment problem is then reduced to the 
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linear programming of minimizing 


n n 


i=l j=l 


( 1 ) 


over all stochastic n x n matrices P := {pi}, i.e. these matrices which satisfy the 2n + linear 
constraints 

n n 

= ; Pi>0 Vi,j G n} . 

k=l k=l 


The Birkhoff Theorem assures us, to our advantage, that the optimal solution of this continuous 
assignment problem is also the solution of the deterministic version. 


The probabilistic version seems to be more difficult since it involves a search on a much larger set 
of n X n stochastic matrices. On the other hand, it has a clear advantage since it is, in fact, a linear 
programming which can be handled effectively by well developed algorithms for such problems. 


In many cases the probabilistic version cannot be reduced to the deterministic problem. For 
example, if the number of sources n and number of targets m not necessarily equal, or when not 
all sources must find target, and/or not all targets must be met, then the constraints are relaxed 
into Y27=iPi — ^ cind/or shall not deal with these extension in the current paper, 

except, to some extent, in section below. 


1.1 Prom the discrete assignment problem to the continuum OMT 

Let p he a probability measure on some measure space X, and v another probability measure on 
(possibly different) measure space Y. Let c = c(x, y) be the cost of transporting x to y. The 
object of the Monge problem is to find a measurable mapping T : X ^ Y which generalizes the 
deterministic assignment perturbation r described above in the following sense: 

T^p = u namely p(T~^{B)) = ^{B) (2) 

for every measurable set B <ZY. The optimal Monge mapping (if exists) realizes the infimum 

inf / c{x,T{x))p{dx) . 

Jx 

The relaxation of Monge problem into Kantorovich problem is analogues to the relaxation of the 
deterministic assignment problem to the probabilistic one: Find the minimizer 

c{p,v) := min / / c{x,y)7r{dxdy) (3) 

■7ren))(^,y) Jx Jy 

among all probability measures tt G IlJ^{p, v) := 

{ Probability measures oa X xY whose X (resp. Y) marginals are p (resp. u)} . (4) 

In fact, Kantorovich problem is just an infinite dimensional linear programming over the huge 
set n^(//, u). The Monge problem can be viewed as a restriction of the Kantorovich problem to the 
class of deterministic probability measures in I\J^{p,v), given by TT{dxdy) = p{dx)5y_T(x) where 
T^/i = u. It turns out, somewhat surprisingly, that the value c{p, u) of the Kantorovich problem 
equals to the infimum Q of Monge problem, provided c is a continuous function on X x K and p 
does not contain a Dirac 5 singularity (an atom) [^. 
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1.2 Semi-finite approximation- The middle way 

Suppose the transportation cost c = c{x, y) on X x Y can be obtained by interpolation of pair of 
functions on X x Z and on Z x T, where Z is a third domain and the interpolation means 

c(x, y) := inf (x, z) -|- (z, y) . (5) 

zSZ 

A canonical example for X = Y = is c{x,y) = c{\x — y\) where c{w) = p > 1. Then ^ is 
valid for Z = and both = 2P~^\w\p. So 

c{x, y) := \x — yY’ = inf \x — z\^ + \z — y\^ ( 6 ) 

zeM** 

for any x,y £ provided p > 1. Note in particular that the minimizer above is unique, z = 
(x + y)/2, provided p > 1, while z = tx + {1 — t)y for any t £ [0,1] if p = 1 . 

Let Z = Zm := {zi,... Zm} C Z is a finite set. Denote 

:= min z) + y) > c(x, y) (7) 

Z&Zm 

the (Zm) semi-finite approximation of c given by (§. 

An optimal transport plan for a semi-discrete cost 0 is obtained as a pair of m—partitions 
of the spaces X and Y. An m—partition is a decomposition of the the space into m mesnrable, 
mutnally disjoint snbset. It turns out that c'^™(/x, z/) can be obtained as 

c^'^{p,iy)= inf / c^^\x, z)p{dx) -£ / c^‘^\z,y)v{dz) ( 8 ) 

{A4,{B4 ^ Ja, Jb, 

ZkzZjm 

where the infimum is on the pair of partitions {Az} of X and {Bz} of Y satisfying p{Az) = v{Bz) 
for any 2 : G Zm- The optimal plan is, then, reduced to m plans transporting Az <£ X to Bz <£Y, 
for any 2; £ Zm, where {Az,Bz} is the optimal partition realizing Q. 

The real advantage of the semi-discrete method described above is that it has a dual formulation 
which convert the optimization Q to a convex optimization on M™'. Indeed, we prove that for a 
given Zm C Z there exists a concave fnnction rfi, 7 : MX —)• M such that 

max 7 (p) = (n, ly) 

and, under some conditions on either p or n, the maximizer is unique up to a uniform translation p -£■ 
p + /3(1,... 1) on M”*. Moreover, the maximizers of 7 ^ yield the unique partitions {Az, Bz; z £ 
Zm} of 0 . 

The accuracy of the approximation of c{x,y) by c^^{x,y) depends, of course, on the choice of 
the set Zm- In the special (but interesting) case X = Y = Z = and c{x,y) = |x — y|'^, a > 1 
it can be shown that c^"*(x, y) — c{x,y) = for any x,y in a compact set, where Zm are 

distributed on a regular grid containing this set. 

From 0 and the above reasoning we obtain in particular 

c^^{p,iy) - c{p,u) >0 (9) 

for any pair of probability measures, and that, for a reasonable choice of Zm, ([^ is of order 
if the supports of y, v are contained in a compact set. 
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For a given m G N and pair of probability measures /r, v and , the optimal choice of is the 
one which minimizes Q. Let 

:= inf - c{n,v) > 0 (10) 

where the infimum is over all sets of m points in Z. Note that the optimal choice now depends 
on the measures themselves (and not only on their supports). A natural question is then to 
evaluate the assymptotic limits 

(^(/r, v) := lim sup v) ; v) := lim inf m^/‘^c™'(|U, u) . 

m^oo m—>-oo 

Some preliminary results regarding these limits are discussed in this paper. 

1.3 Numerical method 

The numerical calculation of Q we advertise in this paper apply the semi-discrete approximation 
(.Zm q£ order m. It also involves discretization of into atomic measures of finite support (n). 
The level of approximation is determined by the two parameters: The cardinality of the supports 
of the discretized measures, n, and the cardinality of the semi-finite approximation m of the cost. 
The idea of semi-discrete approximation is to choose n much larger than m. As we shall see, the 
evaluation of the approximate solution involves finding a maximizer to a concave function in m 
variables, where the complexity of calculating this function, and each of its partial derivatives, 
is of order n. A naive gradient descent method then result in 0{m) iterations to approximate 
this maximum, where each iteration is of order mn. This yields a complexity of order 0{m?‘n) to 
obtain a transport plan on the approximation level of This should be compared to the 

complexity of the Hungarian algorithm [^. We shall not, however, pursue a rigorous complexity 
estimate in this paper. 


1.4 Structure of the paper 


In sectionj^we consider optimal partitions in the weak sense of probability measures, as Kantorovich 
relaxation of solutions of the optimal transport in semi-discrete setting. We formulate and prove a 


duality theorem (Theorem 2.1) which yields the relation between the minimizer of the OMT with 
semi-discrete cost to maximizing a dual function H of m variables. 

In section we define strong partitions of the domains, and introduce conditions for the unique¬ 
ness of optimal solution and its representation as the analogue of optimal Monge mapping. The 
main results of this section is given in Theorem 3.1 In section]^ we introduce an interesting appli¬ 
cation of this concept to the theory of pricing of goods in Hedonic markets, and remark on possible 
generalization of optimal partitions to optimal subpartition. This model, related generalizations 
and further analysis will be pursued in a separate publication. 

In section we discuss optimal sampling of fixed number of centers (m). In particular we 
show a monotone sequence of improving semi-discrete approximation by floating the m centers 


into improved positions. In section 5.2 we provide some assymptotic properties of the error of the 
semi-discrete approximation as m 


oo. 


In section we introduces a detailed description of the algorithm on the discrete level. 
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In section we show some numerical experiments of calculating optimal partitions in the case 
of quadratic cost functions on a planar domain. 

The numerical method we propose in this paper has some common features with the approach 
of Merigot [^, see also [^, as we recently discovered. We shall discuss this issues in section]^ 

1.5 Notations and standing assumptions 

1. X,Y are Polish (complete, separable) metric spaces. 

2. A1+(X) is the cone of non-negative Borel measures on X (resp. for Y). 

3. The weak — * topology on AI+(X) is the dual of Cb{X), the space of bounded continuous 
functions on X (resp. for Y). 

4. Mi{X) is the cone of probability (normalized) non-negative Borel measures in A4+(X) (resp. 
for Y). 

5. For // e Mi{X), u e Mi{Y), 

n](^(/x, u) := {vr G A4i(X x Y) ; /x is the X marginal and v is the Y marginal of tt} 

6. The m-simplex Em := {s := (si,...Sm), Si > 0, = 1} ^ 


2 Optimal partitions 

Definition 2.1. 


i) A m—partition of a pair of a probability measure fi G Afi(X) subjected to f G Em is given by 
m nonnegative measures pLz G Af+(X) on X such that Ylz&Zm h-z = h o-nd dfiz = Vz- The 
set of all such partitions fl := (^i,... pim) is denoted by Vx{p). 

a) If, in addition, v G A4i(T) then {fl,il) G V^{p,,v) iff ft G 'Px{pf) and v G 'Py{v) for some 
r G Em- 


The following Lemma is a result of compactness of probability Borel measure on a compact 
space (see e.g. [^). 

Lemma 2.1. For any r G Em, the set of partitions is compact with respect to the (C*)™’(X) 
topology. In addition, 'P^(/x, z/) is compact with respect to {C*)^{X) x {C*)"^{Y) topology. 

Lemma 2.2. 

C^’"(/X,zx)= min / c^^\x, z)fizidx) + / c^‘^\z,y)vz{dy 

fY Ux Jy 

z^Zjrn 

where c^^ {pi, u) as defined by i0 and {ft, u) G {h, z^) • 


Proof. First note that the existence of minimizer is obtained by Lemma 2.1 
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Define, for z G 

Tz := {{x,y) e X xY; c^^^(x,z) + c^'^\z,y) < c^^{x,y)} C X xY 

such that is measurable in X x D, n = 0 if 2 : 7 ^ and YlzeZm = X x Y. Note that, 
in general, the choice of {F^} is not unique. 

Given vr G FI^(/x, u), let vr^ be the restriction of vr to F,^. In particular Yhz&Zm ~ 
the X marginal of tt^ and Vz the y marginal of Then (/T, u) defined in this way is in V\{yL, v). 
Since by definition c^^{x,y) = c^^\x,z) + c^^\z,y) a.s. vr^. 



X JY 


c^'^{x,y)TT{dxdy) = ^ / / c^™-{x,y)T:z{dxdy) 


Z^Zrr 


X JY 


= ^ / / {c^^\x,z)n:,{dxdy) + / {c^‘^\z,y)nz{dxdy)-Kz{dxdy) 

'y J X J Y J X 

Zkz.Zj'm 

= ^ f c‘^^\x,z)nz{dx) + [ c^‘^\z,y)v^{dy) 

Z&Zm 

Choosing tt above to be the optimal transport plan we get the inequality 

Y'^{y,v)> inf Y] / c^'^\x,z)nz{dx)+ c^^\z,y)vz{dy 


( 11 ) 


To obtain the opposite inequality, let (/U, F) G 'P^(/i, z^) and set := J^dyz = fydi^z- Define 
Tr{dxdy) = YzeZm '''z^f^z{dx)i'z{dy). Then vr G n^(|U,z/) and, from (^7^ 



X JY 


Z^Zrr 


Y'^{x,y)TT{dxdy) = Y^{x,y)rY lJz{dx)vz{dy) 


X Jy 


- [ ic^^\x,z)+ c^‘^\z,y))r^^Hz{dx)uz{dy) 

zeZm 

^ f c‘^^\x,z)nz{dx)+ f c^‘^\z,y)vz{dy) 

r? X Jy 

and we get the second inequality. 


Z^Zm 


( 12 ) 

□ 


Given p = {pz^,. ..pz^)eW^, let 

^zliP’^) ■= ““ c^^Hx,z) +pz ; ^zliP^y) ■= +Pz (13) 

“■ Z&Zm Z&Zm 

=M™(^ •= / &JP^x)p{dx) ■ E^^{p):= f (^zliP^yMdy) . (14) 

J X J Y 

(p) := (p) + Hf- i-p) . (15) 

Lemma 2.3. // // G Ali(X) then for any r G Dm, 

(-H^™)*(-F) := sup (p) - p ■ f = ip, Y VzSz] = min Y] / Y{x,z)pz{dx) . 

\ zeZm / Z&Zrr, ^ 

(16) 
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Analogously, for v G M.i{Y) 


(—:= sup Hf™ (^ — p • r = 

pSK"* 



■ (17) 


Herep-f:= Y.z&z,^^zPz- 


Proof. This is a special case of the general duality theorem of Monge-Kantorovich. See, for example 


22 . It is also a special case of generalized partitions, see Theorem 3.1 and its proof in 24 


□ 


Theorem 2.1. 


sup 


(18) 


Proof. From Lemma 2.2, Lemma |2.3| and Definition |2.1| we obtain 


c^^{p,u)=M [{-E^^^y{-r) + {-E^^y{-r)] . 


(19) 


Note that (—H^"*)*, (—Hf™)* as defined in ( 16, 17), are, in fact, the Legendre transforms of — 

—respectively. As such, they are defined formally on the whole domain M™' (considered as the 
dual of itself under the canonical inner product). It follows that (—H^’")*(r) = {—E^^y{f) = oo 


for r G — Em. Note that this definition is consistent with the right hand side of ( ^ 17), since 
'Pxil^) = 7^y(^) = 0 for f 0 Em- 


On the other hand, and are both finite and continuous on the whole of 
Fenchel-Rockafellar duality theorem (see [22| - Thm 1.9) then implies 

sup 


+ E^^{-p) = jnf (-Hf"‘)*(r) + {-E^’-y{r) . 


rSM" 


The 


( 20 ) 


The proof follows from (15, 19). 
An alternative proof: 


Wecan prove (18) directly by constrained minimization, as follows: (/7, F) G Vy{p,v) iff F(p, 0, V’) := 


Pi(j dpi- j duij + j (fix) f pidx) - ^ Pzidx) j + j ifiy) f uidy) - ^ fzidy) j < 0 

Z^Zjn ^ ^ \ Z^ZjYi / ^ \ Z^Zm / 

for any choice of p G M”^, (f G C’(A), y G C{Y). Moreover, sup^,^ ,^ F = oo unless (/x, F) G V^ip, v). 
We can then obtain from Lemma 2.2 c^™(/x, v) = 


inf sup 

{u^&M+{X),Vz&M+(y)} ,(l,^C(X),'ip&C(Y) 


^ / c^^\x,z)pzidx)+ f c^‘^\z,y)vydy) 

z&Zm 


+Fip,(t),y) 


^ [ (^c^^'>ix,z)+p^-(f)ix)^ p^dx) 

'■&Zm 

+ ^ ^ (c^‘^\z, y)-pz- y{y)) Mdy) + (fpidx) + yv{dy) . (21) 


sup inf 

peiR’",<AeC(x),V'eC(y) {u^&M+(x),i'z&M+(X)} 
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We now observe that the infimum on {fj,z,^z} above is —oo unless z) + Vz — 4>ix) > 0 and 

(z, y) +pz — 'tpiy) > 0 for any z G Zm- Hence, the two sums on the right of (21) are non-negative, 


so the inhmum with respect to ’^z} is zero. To obtain the supremum on the last two integrals 
on the right of © we choose as large as possible under this constraint, namely 


(b(x) = min z) -|- Pz 


■4){y) = min c^‘^\z,y) - pz 

Z&Zm 


SO 4>{x) = V’(y) = Czii-P^y) by definition via (|l3|). 


_ f( 2 ) 


□ 


3 Strong partitions 

We now define strong partitions as a special case of partitions (Definition |2.1[ ). 

Definition 3.1. 

i) A partition p, G V^ip) is called a strong m—partition if there exists m measurable sets Az C X, 
z G Zm which are essentially disjoint, namely p{Az n Ap) = 0 for z z and p{Uz^Zm^z) — 
X, such that pz is the restriction of p to Az. The set of strong m—partition corresponding to 
r G Sm is denoted by Vx{p). 

a) In addition, for v G AAiiY) then {p,u) G V^ip^i^) iff P ^ ^ ^ ^y(^) some 

r G Sm- In particular, a strong m—partition is composed of m p measurable sets Az (Z X 
and m v measurable sets Bz CY such that dp = dv for z G Zm- 

Assumption 3.1. . 

a) p € All (A) is atomless and p{x-,c^^\x,z) — c^i^\x,z) = p) = 0 for any p G M and any 

Z,Z G Zm- 

b) u £ Afi(y) is atomless and v{y]c^‘^\z,y) — c^‘^\z',y) = p) = 0 for any p G M and any 

Z,Z G Zm- 


Let us also define, for p G M”* 

Az{^:={x£X-, c(b(x,z)-Kp^ = 4^j^(p,x)} ; Bz{p) := {y ^Y; c(2)(z,y)-Fp^ = 4^(p,p)} . 

( 22 ) 

Note that, by (13, 14) 

(23) 

(24) 


likewise 




-^Zr. 


Y] [ Y^\x,z) +Pz)pidx) 

ip)=Y I +Pz)’^idy) ■ 

JBPp) 


Lemma 3.1. Under assumption 3.1 (a) (resp. (b)) 


i) For any p G M™, {Az{p)} (resp. {Bz{p)}) induces essentially disjoint partitions of X (resp. 








a) 


(resp. is continually differentiable funetions on 


Q’^Zm az:Zm 

^ = p{Affp)) resp. = v{Bffp)) . 


dpz dp. 

This Lemma is a special case of Lemma 4.3 in [W]. 


Theorem 3.1. Under either assumption 3.1-(a) or (b) there exists a unique minimizer ro of (19). 
In addition, there exists a maximizer pQ G M™' and either (in case (a)) { 742 (^ 0 )} or (in 

ease (b)) {Bz{—po)} induces a eorresponding strong m—partition in (a) 'P^(/x) or (b) Vyiv). In 
particular, i/both (a+b) holds then {^ 2 (^ 0 ); {-^^(—po)} induees a strong m—partition inV^ih^^)} 
and 

TToidxdy) := ^ {ro^z)~^lA,(fo){x)lB,(-po){y)h{dx)i'{dy) (25) 

zeZm-,ro,z=ui^ziPo)) 

is the unique optimal transport plan for u). 


Proof. Note that '^{p) — r-p is invariant under additive shift for H = and r G Indeed, 

H(p + al) = + a for any a G M where 1 := (1,... 1). So, we restrict the domain of H to 

pG i?™ , p - 1 = 0 . (26) 

Assume (a). Given r G Assume first 

G (0,1) for any z e Zm . (27) 

We prove the existence of a maximizer pQ, 

(-S^™)*(-^) = S^™(Po) -Po • r > r 

for any p G M™'. Let Pn be a maximizing sequence, that is 

lim 'Efff{pn)-Pn • r = {-r.fff)*{-r) 

(c.f. (0). 

Let ||^|2 := (E zeZmPz'^^^'^ be the Euclidian norm of p = {pzi, ■ ■ - Pz^) ^ If prove 

that for any maximizing sequence Pn the norms ||pn ||2 are uniformly bounded, then there exists a 
converging subsequence whose limit is the maximizer po- This follows, in particular, since is 
a closed (upper-semi-continuous) function. 

Assume there exists a subsequence along which ||pn ||2 —^ 00. Let pn '■ = Pn/\\Pn\\2- Let 


'^friPn) -Pn-r-= [^fr{Pn) “ Pn ' (pn)] + Pn ' (Vj^“^-(pn) - f) 

~ ~ Pn ■ {Pn)^ + ||Pn||2Pn ‘ (Pn) ~ ?’) • (28) 

In addition, by (23) and Lemms[3d^(ii) 


oo < / vam. cd\x, z)p{dx) < {p) — p ■ V{f)j\ = / cd\x,z)p{dx) 

JX Ja^p) 

< / max c*'^^(x, 2 ;)/r(dj;) < OO . ( 29 ) 

Jx Z&Zm, 
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(30) 


By (28- 29) we obtain, for 


oo, 


lim pn 
n^oo 


- f) = 0 


Since pn lives in the unit sphere S'™ ^ in M™ (which is a compact set), there exists a subsequence 
for which pn po := (po,zi, • ■■Po,z^) G S™“^ Let P_ := and J_ := {z G Zm ]Po,z = 

P_}. 

Note that for n —)■ oo along such a subsequence, pn^z — Pn,zr —^ —oo for z G J-, z/ ^ J-. It 
follows that Az,{pn) = ^ if zf ^ J- for n large enough, hence U^gj_A^(p„) = X for n large enough. 
Let be the restriction of p to Az{pn)- Then the limit p"^ pz exists (along a subsequence) 
where n —)• oo. In particular, by Lemma |3.1| 


lim 

n^oo 


dPn,z 


(Pn) 



dpz 


while Pz ^ 0 if only if z G J-, and Ylz&j^ Pz = P- Since po^z = P- for z G J- is the minimal value 
of the coordinates of po , it follows that 


lim Pn ■ {V^E^^ipn) -f) = -r-po + Y] / dp, 

n^oc V i' / I 

z&J- 


= -r- Po + P- ■ 


Now, by (27), r • po > P- unless J_ = Zm- In the last case we obtain a contradiction of (26) 
since it implies po = 0 which contradicts po G S™“^. If J_ is a proper subset of Zm we obtain a 
contradiction to (30). 


If (27) is violated we may restrict to domain of to a subspace by eliminating all coordinates 
z G Zm for which = 0. On the restricted subspace we have a minimizer pQ by the above proof. 
Then we may exte nd po by assigning pz sufficiently small if = 0. This guarantees Az{po) = 0, 
hence (Lemma 3.1^ dE^"^ /dpz = 0 for any such Hence the extended po is still a critical point of 


— f- p, and is a maximizer by concavity of . 


Next, we prove that Az{po) is a unique optimal partition of X. Let p G be a minimizer of 


(16). Since dpz = Vz, Pz = P, (16) implies 


Y [ c^^Hx,z)pz{dx) = {-EY*i-'^ 


and 


“^'^0 = / CzliPo^x)dp-Po-r= Y [ {^zliPo^^)-Po,z) dpz , 

Z&Zrr, 


SO 


X] / -po^z-c^^\x,z)^ Pz{dx) = f) . 

Z&Zrr, 


On the other hand, ^^^^{pq,x) — po.z — c^^\x,z) < 0 for any x G X hy definition (13), so we must 
have the equality 


^zliPo^x) =po^z + c^^Kx,z) 


(1)/ 


10 













a.e. on supp(//^). Hence supp(/i^) C Az{po). Since Az{po) are mutually disjoint and YlzeZm ~ 
then pz is necessarily the restriction of p to Az{po). On the other hand, for any p ^ po mod Ml 
there exists z G Zm for which p (^Az{po)AAz{^'^ / 0. This implies that the strong partition A{po) 
is the unique one. 


The same result is applied to (p) — Po ' If we show that the minimizer Tq of the right 
side of (20) is unique, then it follows that the maximizer pQ of the left side of (20) is unique as well 
(up to IM), and, in particular, the optimal partition is unique. Hence, we only have to show the 
uniqueness of the minimizer of the right side of (20). This, in turn, follows if either (— 

V 


or 


is strictly convex. 

To prove this we recall some basic elements form convexity theory (see, e.g. [BC]): 


i) If T is a convex function on M™ (say), then the sub gradient dF at point p G M™ is defined as 
follows; q G dF{p) if and only if 

F (pt) — F{p) > q ■ {pf — p) \/pf G M™' . 


ii) The Legendre transform of F: 

F*{^ := sup p-q-F{p) , 

pSK"* 

and Dom{F*) C M™ is the set on which F* < oo. 

iii) The function F* is convex (and closed), but Dom{F*) can be a proper subset of M™ (or even 

an empty set). 

iv) The subgradient of a convex function is non-empty (and convex) at any point in the proper 

domain of this function (i.e. at any point in which the function takes a value in M). 

v) Young’s inequality 

F{p) + F*{q) >p-q 

holds for any pair of points (p, q) G M”^ x M™'. The equality holds iff g G dF{p), iff p G dF*{q). 

vi) The Legendre transform is involuting, i.e F** = F if F is convex and closed. 

vii) A convex function is continuously differentiable in the interior of its proper domain iff its 
subgradient at any point in the interior of its domain is a singleton. 


Returning to our case, let F := It is a closed, convex, proper and continuously differentiable 

function defined everywhere on M™. Assume (—is not strictly convex. It means there exists 
fij^r 2 & Dom{—E^^)* for which 




ri+r2, i-^n^Tin) + {-^n^)*{r2) 


Let f := fi/2 + f2/2, and p G Then, by (iv, v) 


(31) 


0 = (-H^™)*(iO + {-E^-y*{p)-p- f= (-H^- 


)*(r)- 


-p-r . 


(32) 
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By (31 32): 


J fi) + ^ ((-H^'")*(r2) - ^fr{^-p- fs) = 0 


while (v) also guarantees 


It follows 


{-'^l"")*{ri)-'^l'"{p)-p-ri>Q , i = l,2. 

-^n'"{p)-p-ri = 0 , i = l,2, 


SO, by (v) again, {ri,r 2 } E dE^^{p). This is a contradiction of (vii) 


since 




is continuously 


differentiable everywhere on by Lemma 3.1 


Finally, we prove that ttq given by (25) is an optimal plan. First observe that ttq G Il{p,u), 
hence 

/ / c^’"{x,y)7ro{dxdy) . 

Jx Jy 

Then we get, from ([^ 

c^™(/r,j/)< / / c^^{x,y)7ro{dxdy) < V / {c^^\x, z)p{dx) + c^ 2 ){y, z)i'{dy)) 

JxJy ^^z^Ja,(po)xB,(-po} 

= V ( / c^^\x,z)p{dx) + [ c^^\z,y)i^{dy)] =E{po) < 

zeZm V^-^Po) dB4-po) J 

where the last equality from Theorem 2.1, In particular, the first inequality is an equality so vro is 
an optimal plan indeed. □ 


4 Pricing in hedonic market 


In adaptation to the model of Hedonic market there are 3 components: The space of consumers 
(say, X), space of producers (say Y) and space of commodities, which we take here to be a finite 
set Zjn := {zi,... Zm}- The function := z) is the negative of the utility of commodity 

z G Zm to consumer x, while := c^^\z,y) is the cost of producing commodity 2; G Zm by the 
producer y. 


Let // be a probability measure on X representing the distribution of consumers, and n a 
probability measure on Y representing the distribution of the producers. Following we add the 
’’null commodity” zq and assign the zero utility and cost c^^\x,zo) = c^‘^\zo,y) = 0 on X (resp. 
Y). We understand the meaning that a consumer (producer) chooses the null commodity is that 
he/she avoids consuming (producing) any item from Zm- 


The object of pricing in Hedonic market is to find equilibrium prices for the commodities which 
will balance supply and demand: Given a price pz for z, the consumer at x will buy the commodity 
z which minimize its loss c^^\x,z) + Pz, or will buy nothing (i.e. ’’buy” the null commodity zq) 
if min^g^^ c(^)(x, z) + Pz > 0), while producer at y will prefer to produce commodity z which 
maximize its profit —c^‘^\z,y) +Pz, or will produce nothing if maXz^Zm ~c^‘^\z,y) +pz < 0. Using 


notation (13 15) we define 


a;) := min{4’^(p, x), 0} ; ^^{p, y) := min{^f2(^’ 2/)’ 0} 


(33) 
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(34) 


iP) ■= / ixiP^x)ix{dx) 


lx 


-^^p) :=^. 


iY{PiV)^{dy) 




i-p • 


(35) 


Thus, is the difference between the total loss of all consumers and the total profit of 

all producers, given the prices vector p. It follows that an equilibrium price vector balancing 
supply and demand is the one which (somewhat counter-intuitively) maximizes this difference. The 
corresponding optimal strong m—partition represent the matching between producers of {Bz C T) 
to consumers {Az G X) oi z € Z. The introduction of null commodity allows the possibility that 
only part of the consumer (producers) communities actually consume (produce), that is Uz^Zm^z C 
X and UzeZmBz C Y, with Aq = X - Uzez^Az (Bq = Y - UzeZm^z) being the set of non-buyers 
(non-producers). 


From the dual point of view, an adaptation Cq"* (x, y) := min{c^"* (x, y), 0} of Q (in the presence 
of null commodity) is the cost of direct matching between producer y and consumer x. The optimal 
matching (Az,Bz) is the one which minimizes the total cost over all su5-m—partitions 

V"^{pi,v) as defined in Definition 


3.1 


(ii) with the possible inequality p{GAz) = v{GBz) < 1. 


5 Dependence on the sampling set 


So far we took the smapling set C to be fixed. Here we consider the effect of optimizing Z^ 
within the sets of cardinality m in Z. 

As we already know from ([^[^, c^^{x, y) > c{x, y) on X xY for any {x,y) ^ X xY and Z^ C Z. 
Hence also c'^™(/x, u) > c{p, v) for any y, G AIi and any Z^ C Z as well. An improvement of Z^ 
is a new choice Z^'^ C Z of the same cardinality m such that {p, v) < v). 


C Z, once the optimal partition is 
the improvement depends on the measure p, v. 


In section 5.1 we propose a way to improve a given Z, 
calculated. Of course 


In section 5.2 we discuss the limit m —)• oo and prove some assymptotic estimates. 


5.1 Monotone improvement 


Proposition 5.1. Define SJ) on M™' with respect to Zm '■= {zi,... Zm} ^ Z as in (15). Let 
{p,il) G Vx{p,iz) be the optimal partition corresponding to c^'^{p,u). Let ({i) ^ Z be a minimizer 
of 

ZbC^ f c^^\x,C)hzi{dx) + [ c‘'^'>{C,y)i'zi{dy) . (36) 

Jx Jy 

Let Zffi^ := {C(l), • • • C("^)}- Then c^^'" {p, v) < c^"‘{p, v). 


Corollary 5.1. Let Assumption 3.1 (a+b), andpo be the minimizer of inMX. Let {Az{po), Bz{—po)'\ 
be the strong partition corresponding to Zm as in (22). Then the components of are obtained 
as the minimizers of 


Z 3 ( i-j. c^^\x,C)p{dx) + 

lYiPo) 


Bz(-po) 


P‘^HC,y)i^idy) 
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2.1 




{p) < := maxjRm for any p G 


proposition 5.1): Let be defined with respect to Z^'^. By Lemma 2.2 and 


Proof, (o: 

Theorem 

so inaxKm r.i_l [p) = [p, v) < max]R™ [p) = v). 

Remark 5.1. If c is a quadratic cost then is the center of mass of Az{po) and Bz{—po): 

K^ziPo)) + i'{B^{-po)) 

We shall take advantage of this in section \Kl\ 

Let 

c^(p,u) := inf c^’"(/x, z/) . 

ZmCZ ; #iZm)=m 


□ 


Let Z!f^ := {z^,... C .Z’ be a sequence of sets such that is obtained from Z^ via (36). 
Then by Proposition |5.1| 


ryk 

c{p, v) <d^{p,v) < ■ • • c”" {p, f) < {p,v) < ... {p, v) . 

Open problem: Under which additional conditions one may gurantee 

lim c^^(/r, z^) = c™'(/r, z/) ? 

fc^OO 


5.2 Assymptotic estimates 


Recall the definition 


( 10 ) 


(Tibw) ■■= inf c^"*(ir,z^) -c(/r,z/) > 0 

ZimC.Zi 

Consider the case X = Y = Z = W^ and 


c(x, y) = min h{\x — z\) + h{\y — z\) 

where h : —)■ M'*' is convex, monotone increasing, twice continuous differentiable. Note that 

c{x,y) = 2h{\x -y\/2). 

Lemma 5.1. Suppose both p and v are supported on in a compact set in M'^. Then there exists 
C = C{p, i^) < oo such that 

hmsupm^/‘^0™(/z, z/) < z/) . (37) 

m^oo 


Proof. By Taylor expansion of z —>■ h{\x — z\) + h{\y — z\) at zq = {x + y)/2 we get 

h{\x - z|) + h{\y - z\) = 2h{\x - y\/2) + - v) ' - ^o)]^ + o^{z - zo) . 

Let now Z^ be a regular grid of m points which contains the support K. The distance between 
any z ^ K to the nearest point in the grid does not exceed for some constant C{K). 
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Hence Cmix,y) — c{x,y) < sup \h''\C{K)'^m if 2;, 2/ £ K. Let 7rQ{dxdy) be the optimal plan 
corresponding to y, v and c. Then, by definition, 


c{y,v)= / / c{x,y)'Ko{dxdy) ; < / / Cm{x,y)'Ko{dxdy) 

Jx Jy Jx Jy 


4>^{y,iy)< / / {cm{x,y) - c{x,y))'KQ{dxdy) <sup\}i'\C{Kf‘m 

Jx Jy 

since tto is a probability measure. 




□ 


If h[s) = 2°" ^5°^ (hence c{x,y) = \x — y\^) then the condition of Lemma 5.1 holds if cj > 2. 


Note that if /r = then c{y, y) = 0 so 4>'^{y,y) = iniz^^z y)- In that particular case we 

can improve the result of Lemma |5.1| as follows: 


Proposition 5.2. If c{x,y) = \x — y\^, (t> 1, X = Y = Z = W^ and u = y = f{x)dx 


lim = Cd,a 

m^oo 


{d-\-G)/d 


(38) 


where Cd,a is some universal constant. 


Proof. From (15), (p ) 

be p = 0. By Theorem |2.1| 


Using (13, 14) with c^^'>{x,y) 


= H^"*(p) + is an even function. Hence its maximizer must 

= c^‘^\y, x) = 2^~^\x — y\^ we get 

= 2^ I min |a: — z\^yL{dx) . 

J^dZGZm 


We then obtain (38) from Zador’s Theorem 10, 25, 10 


□ 


Note that Proposition 5.2 does not contradict Lemma |5.1[ In fact cr > 2 it is compatible with 


the Lemma, and (37) holds with C{y,y) = 0 if u > 2. If u G [1, 2), however, then the condition of 
the Lemma is not satisfied (as h is not bounded near 0), and the Proposition is a genuine extension 
of the Lemma, in the particular case y = v. 

We can obtain a somewhat sharper result for any pair y, iz in the case <7 = 2, which is presented 
below. 


Let X = Y = Z = c{x,y) = \x — yp, /x, z/ are Borel probability measures which admits a 
finite second moment. Assume /x is asbolutely continuous with respect to Lebesgue measure on W^. 
In that case, Brenier Polar factorization Theorem implies the existence of a unique solution to 
the quadratic Monge problem, i.e a Borel mapping T such that T^/x = iz. Let A = f{x)dx be the 
McCann interpolation between y and iz, that is, A = {I/2 + T/2)^y. We know that A is absolutely 
continuous with respect to Lebesgue as well. 


Theorem 5.1. Under the above assumptions, 


limsupm^/'^</>™'(/x, zx) < 4.Cd,2 
m—^oo 




{d+2)/d 
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Proof. Let S be the optimal Monge mapping transporting A to u, i.e. S^X = i/ is a solution of 
Monge problem 

f \S{x) — x\‘^X{dx) = min f \Q{x) — x\^X{dx) . 

Note that if y = {T{x) + x)/2 then S{y) = T{x). Then, since A = (1/2 + T/2)^y, 


\Siy)-yfXidy) = 


T{x) — X 


fi{dx) = c{n, v)/^ . 


Also, if y = {T{x) +x )/2 then 2y — S{y) = x. If follows that 2/ — S' is the optimal Monge mapping 
transporting A to /x, that is, 

f \S{x) — 2x\‘^X{dx) = f \S{x) — x\‘^\{dx) = min f \Q{x) — x\‘^X{dx) 

jRd J^d Q;Q#A=/ijRd 


SO 


c{y,u) = 2 / |S(x) — xpA(dx) + 2 / \S{x) — 2x\‘^X{dx) = A / |S(x) — xpA(dx) . (39) 

jR<i JRd jRd 


Given z G Zm, let 


14 := G \x — z\ <\x — z \ Vz G 


Since = 1^“ and A(I4 H F/) = 0 for z 7 ^ z then (39) implies 


= / \S{x) - x\‘^X{dx) . 


Let Vz '■= S^A[I4, := {21 — S)^A[I4. Form Lemma 2.2 


‘(/X, z/)<2( ^ j \x - z\^ yz{dx) + \x - zfuzidx)] 

\ZGZm ZGZm / 

= 2 f {|S(a:) — z|^ + |2x — S(x) — zp} A(dx) 

Jvz 


(40) 


(41) 


(42) 


By the identity 


4|z — x\^ = 2 {|S(x) — z|^ + \2x — S{x) — z|^} — 4|S(x) — x\^ 
This, together with (41, 42) and implies 

(/>™'(/x, ix) < 4 ^ f |x — z|^A(fix) . 


(43) 


By (40), fvz “ z\‘^X{dx) = J^d min^g^^ |a: - z\‘^X{dx) := 4>{X,Zm). Since (43) is valid for 

any we get the result from Zador’s Theorem (10, 25, |10|. □ 
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6 Description of the Algorithm 


We now spell out the proposed algorithm for approximating of the optimal plan c(/i, u). We assume 
that c is given by ([^. We hx a large numbers ni,n 2 (not necessarily equal) which characterizes 
the fine sampling, and much smaller m characterizing the partition order. Then we choose an 
appropriate sampling: In X we set Hm ■= YllLi SiSxi for // and on Y we set r'ri 2 := for 


At the first stage we choose := {zi,... z^} G Z™, and dehne 


ni 


n2 


’^oip) ■='^Si mm [c^^\xi,z^A+pj] + '^Ti mm [c^‘^\z'^, ,yi) - pj] 


2=1 


2=1 


l<ji<m 


Next we choose a favorite method to maximize Hq on M"*. It is helpful to observe that Hq is 
differentiable a.e. on M™. Indeed, let 


A^(p) := {i G (l,...ni), z°) + pj = min [c^^\xz, + pk]} 

■’ l<k<m 

Bjip) ■={i n 2 ), c^‘^\z^,y^) +pj = rnin {z'l,yz) + pk]} ■ 

1 <r hxT 


Then 


d^o 

dpj 


Y. Y 


provided Aj{p) n Afc(p) = 0 and Bj{—p) n Bk{—p) = 0 for any k ^ j. 

Let po be the maximizer of Hq on M™, A° := A°(po), Bj := Bj{—po). 

At the I step we are given := {z[,... G Z"^, pi the maximizer of 

ni 712 

^i{p) ■=Y.^^ +Pj] + y^T^ min [c^'^\z],yz) - Pj] 


and the corresponding A^- := A^-{pi), Bj := Bj{—pi) where 

A'-j{p) := {i e ni), Zj) + pj = min [c‘'^\xz, z^) + pk]} 

•’ l<k<m 

Bjip) := {i G (1,. . . 712 ), c^‘^\z^j,yz) +Pj = min [c^‘^\zi,yz) + Pk]} ■ 

\ < k>< -TV) 


We define z^^^ as the minimizer of 


Y SzC^^Kxz,z) + ^ TzC^‘^\z,yz) 

ieAi ieB’ij 


(44) 


and set := G Z^. Now 


rii 


n2 


'~l+i{p) ■.= '^Sz min [c(^)(x;„ 2 :'+^)+pj] + min ,yz) - Pj] ■ 

• ^ 1 ■<« 7 \ 771 ‘ ^ 1 <C 'J <r m 


2 = 1 


2=1 




From these we evaluate the maximizer p/+i the maximizer of and the sets B^A^. 
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Remark 6.1. The maximizer pi at the I stage can be used as an initial guess for calculating the 
maximizer at the next stage. This can save a lot of iterations where the stages where changes 
of the centers —>• is small. 

Using Proposition |5.1| and Corollary |5.1| we obtain a monotone non increasing sequence 
So(po) > • • • > ^i{pi) > Ei^i{pi+i) ...> c(/r, u) . 

The iterations stop when this sequence saturate, according to a pre-determined criterion. 


6.1 Application for quadratic cost 


As a demonstration, let us consider the special (but interesting) case of quadratic cost function 
c{x,y) = \x — y\^ on Euclidean space X = Y. We observe the trivial inequality \x — yp = 
2m\nz&x\\x — zj^ + \y — -zp]- Hence we may approximate \x — y\^ by 

c^’"(x, y) := 2 min [|x — + ly — z|^] > |x — yp . (45) 


So, we use c^^\x, z) := 2|x — Zz\‘^, c^‘^\zz,y) ■= 2|y — 

The updating (44) takes now a simpler form due to Remark 
the center of mass 


5.1 


Indeed, is nothing but 


= 


Yhi&x. A 


7 Some experiments with quadratic cost on the plane 


In this section we demonstrate the algorithm for quadratic cost. The pair (y, X) is always considered 
to be uniform Lebesgue measure on the unit square B := {(xi,A 2 ); 0 < xi,X 2 ) < 1}. It is 
sampled by an empiric measure of regular grid composed on 400 points x[*^ = i/20, x^'^^ = j/20, 
y{{i/20, j/20)}) = 1/400, 1 < i,j < 20. The image space {Y,n) is, again, a probability measure 
on the plane which depends on the particular experiment. The number of centers m = 10 and their 
initial choice is arbitrary within the unit square. 


In the first experiments we used a given mapping T := (Ti,T 2 ) : B —)> M^, and defined 
(y, z/) according to T = T{B), v = T^y. In that case the naturel sampling is just (yf\y 2 '^^) = 
(ri(i/20),T2(jy20)), and v[{{y}i\yf)\) = 1/400. 


In all these experiment we used T^ = d^/dxk, k = 1,2, where 4>(xi,X2) = 0.5(xf + x^) + 
A(cos(xi + 2x2) — sin(xi — X2)). FigsT][2 shows the saturated result for different values of A. 

Fig|^?? show pair of partitions on the X square. The right square is the image under (V<h) ^ 
of the partition in the left square. Note that for small values of A the two partitions looks identical. 
This is, in fact, what we expect as long as $ is a convex function. Indeed, the celebrated Brenier’s 
theorem of polar factorization implies just this! For larger values of A, is not convex and we 
see clearly the difference between these two partitions. 


In the second class of experiments we used different domains for Y (e.g. T shaped, I shaped 
and A shaped) which are not induced by a mapping. Fig. ??-?? display the induced partitions 
after saturation for different initial choices of the centers ^ 2 . It demonstrates that the saturated 
partition may depend on the initial choice of the centers. 


18 







Figure 1: Partition for <5, A = 0.2 




Figure 2: Partition for <f>, A = 1.2 



Figure 3: Comparison of partitions A = 0.05 



Figure 4: Comparison of partitions A = 0.1 



Figure 5: Comparison of partitions A = 0.2 
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8 Comparison with other semi discrete algorithms 


Applications of semi-discrete methods for numerical algorithms where introduced in paper by 
Merigot , followed by a paper of Levy [3 . Here we indicate the similar and different aspects of 
our proposed algorithm, compared to 


The starting point of Merigot-Levy algorithm for quadratic cost involves a discretization Vm 
of the target measure v. For optimal plan for transporting ^ is obtained by 

maximizing 




min 

l<i<m 


m 

x-yi\‘^ + Pi]n{dx) - ^ riPi . 

i 


(46) 


This is equivalent to the function we defined (for the special case of quadratic cost) as 


whose maximum over 


is 


{p)-p-f, 

)*(—r) as defined in (16). The optimal partition induced by 
maximizing (46) is refined by taking finer and finer discretization of v with increasing number of 
points m. The multi-grid method is, essentially, using the data of the maximizer p corresponding 
to Vm as an initialization for the m -|- 1 level maximization corresponding to (46). 


In the present paper we take 
function c = c{x,y) via ([^. It is, 
(in the quadratic case), as we can 
step forward we could reduce the 


a different approcah, namely the semi-discretization of the cost 
in fact, equivalent to a two sided discretization analogus to (46) 
observe from (19). However, by carrying the duality method one 
optimization problem to a single one over 


via Theorem 2.1 
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Figure 6: Comparison of partitions A = 0.5 
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